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ABSTRACT 

Context. In analyses of stellar spectra and colours, and for the analysis of integrated light from galaxies, a homogeneous grid of 
model atmospheres of late-type stars and corresponding flux spectra is needed. 

Aims. To construct an extensive grid of spherically symmetric models (supplemented with plane-parallel ones for the highest surface 
gravities) , built on up-to-date atomic and molecular data, and make it available for public use. 
Methods. The most recent version of the MARCS program is used. 

Results. We present a grid of about 10 4 model atmospheres for stars with 2500 K < T e g < 8000 K, -1 < logg = log(GM/R 2 ) < 5 (cgs) 
with various masses and radii, -5 <[Me/H]< +1, with [cc/Fe] = 0.0 and 0.4 and different choices of C and N abundances, including 
"CN-cycled" models with C/N = 4.07 (solar), 1.5 and 0.5, C/O ranging from 0.09 to (normally) 5.0 to also represent stars of spectral 
types R, S and N, and with 1.0 < £ t ^5 km/s. We also list thermodynamic quantities (T, P g , P e , P, partial pressures of molecules, 
etc) and provide them on the World Wide Web, as well as calculated fluxes in approximately 108,000 wavelength points. Underlying 
assumptions in addition to ID stratification (spherical or plane-parallel) include hydrostatic equilibrium, mixing-length convection 
and LTE. A number of general properties of the models are discussed, in particular in relation to the effects of changing abundances, 
of blanketing and of sphericity. We have found and illustrate positive and negative feed-backs between sphericity and molecular 
blanketing. Models are compared with other available grids and excellent agreement is found with plane-parallel models of Castelli 
and Kurucz (if convection is treated consistently) within the overlapping parameter range. Although there are considerable departures 
from the spherically symmetric NextGen models, the agreement with more recent PHOENIX models is gratifying. 
Conclusions. The models of the grid show considerable regularities, but some interesting departures from general patterns occur for 
the coolest models due to the molecular opacities. We have tested a number of approximate "rules of thumb" concerning effects of 
blanketing and sphericity and found them to often be astonishingly accurate. Some interesting new phenomena have been discovered 
and explored, such as the intricate coupling between blanketing and sphericity, and the strong effects of carbon enhancement on 
metal-poor models. We give further details of line absorption data for molecules, as well as details of models and comparisons with 
observations in subsequent papers. 
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giants 



1. Introduction 

Since the first grids of line-blanketed model-atmospheres for 
late-type stars were published (Carbon & Gingerich 1969; 
Querci, Querci & Tsuji 1974; Peytremann 1974; Gustafsson et 
al. 1975; Kurucz 1979; Johnson, Bernat & Krupp 1980), there 
has been a very impressive improvement in underlying data, in 
particular for atomic and molecular absorption. Accurate con- 
tinuous absorption coefficients for a number of heavy elements 
have been calculated within the Opacity Project (Seaton et al. 
1994), and the Iron Project (Bautista 1997). The Opacity Project 
also contributed accurate transition probabilities for a wealth of 
spectral lines for elements where the LS coupling approxima- 
tion was applicable. Thanks to systematic efforts by R. Kurucz 
(see http://kurucz.harvard.edu) line lists with transition proba- 
bilities for millions of lines of heavy elements have been cal- 
culated. Also, the admirable and systematic work by experi- 
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mental physicists (e.g. Blackwell et al. 1989, Nave et al. 1994, 
Hartman et al. 2003, Lawler et al. 2007) has led to the identi- 
fication and measurements of many more metal transitions, e.g. 
from Fe i and Fe n, than existed before. The results of these ef- 
forts have been made easily accessible in data bases, such as the 
VALD (see Piskunov et al. 1995, Stempels, Piskunov & Barklem 
2001, see also http://www.astro.uu.se/~vald). Additional signifi- 
cant progress has been made in the calculation of accurate damp- 
ing constants for strong atomic lines (Anstee & O'Mara 1995, 
Barklem, Piskunov & O'Mara 2000a) as well as for hydrogen 
self broadening (Barklem, Piskunov & O'Mara 2000b). 

Also for molecular absorption, impressive progress has 
been made, starting from practically nothing 30 years ago (see 
J0rgensen 1994b and J0rgensen 2005 for reviews). Very ex- 
tensive line lists have thus been calculated for most diatomic 
and polyatomic molecules that contribute opacity in stellar at- 
mospheres. The lists are based on laboratory measurements of 
wavelengths and gf values with simple theoretical extensions, 
or are the result of more extensive quantum-mechanical ab-initio 
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calculations. Even if the absorption cross-sections are not always 
of satisfactory quality as yet, the existence of these lists makes 
it possible to calculate models for, e.g., M and C stars that are 
realistic enough to compare reasonably well with observations. 

Simultaneously with these improvements of basic physics, 
steps to improve the physical consistency of the models have 
also been taken. The traditional assumptions of plane-parallel 
stratification in homogenous layers, of stationary hydrostatic 
equlibrium, mixing-length convection and LTE have stepwise 
been possible to relax. Grids of spherically symmetric mod- 
els for giants and supergiants were made by Plez, Brett & 
Nordlund (1992) and J0rgensen, Johnson & Nordlund (1992). 
Dynamic pulsating model atmosheres for cool giants have been 
developed by Wood (1979) and Bowen (1988), and these were 
improved with time-dependent dust formation by Fleischer et 
al. (1992) and Hofner & Dorfi (1997), and further developed 
to include frequency-dependent radiative transfer by Hofner 
(1999) and Hofner et al. (2003). Nordlund & Stein and col- 
laborators (Nordlund 1982, Nordlund & Dravins 1990, Stein & 
Nordlund 1998, Asplund et al. 1999, 2000a) developed 3D sim- 
ulations with proper hydrodynamics and radiation fields taken 
into acount for solar-type stars, and Freytag (2001) developed 
full "star-in-a-box" models for supergiants. These simulations 
show a striking agreement with observations of solar granulation 
and spectral line profiles for solar type stars, and clearly demon- 
strate the qualitative difference between traditional ID models 
and reality: while the temperature structure in the upper layers 
of the ID models is determined by radiative cooling and heat- 
ing, the radiative heating in real late-type stars is balanced to a 
significant extent also by expansion cooling of uprising gas. 

A consistent treatment of radiative transfer in models for 
late-type stars without making the assumption of LTE is compli- 
cated, due to the great number of atomic and molecular species 
affecting the radiative field, the wealth of levels and transitions 
in these species, and the lack of basic data for these transitions - 
not the least the lack of cross sections for atomic and molecular 
collisions with electrons and hydrogen atoms. From the first at- 
tempt to construct a reasonably realistic NLTE model for a late- 
type star (the Sun) made by Anderson (1989), the development 
of algorithms and computers has now made it possible to calcu- 
late grids of such models (Hauschildt et al. 2002). The lack of 
accurate collision cross sections, however, is still a major prob- 
lem in these efforts. 

In the present paper we present an extensive grid of model 
atmospheres for late-type stars. The models are still classical, 
in the sense that they are one-dimensional, i.e. spherically sym- 
metric or plane-parallel, with LTE assumed and convection taken 
into account by using the standard mixing-length theory. In spite 
of the inadequacy of these underlying physical assumptions, 
which must be abandoned in accurate work as soon as possible, 
we still think the present grid will be useful for years to come, for 
various applications extending from studies of individual stars to 
galactic evolution and populations in external galaxies. 

The grid extends from models for A-type stars to M, S and C 
star models, from dwarfs to supergiants in luminosity, and from 
10~ 5 times solar to 10 times solar in metallicity, and with various 
choices of parameters like radius and mass, abundances of C, N 
and O, and of "a elements" (Ne, Mg, Si, S, Ar, Ca, Ti) relative 
to Fe, as well as of microturbulence parameter. We also provide 
model Spectral Energy Distributions (SEDs), sampled in about 
10 5 frequency points. As regards methods and underlying data 
used, the model grid is homogeneous. 

The grid is presented in a number of papers. In this initial 
paper we first (Sec. 2) sketch the development of our computer 



code MARCS, from which models produced at different earlier 
stages are in wide use, give a general outline of the methods now 
used (Sec. 3 and 4), and underlying data (Sec. 5), as well as prop- 
erties of the grid (Sec. 6). Models are also compared with models 
from other contemporary grids (Sec. 7). In subsequent papers we 
shall discuss models for A-G stars (Edvardsson et al., Paper II), 
K and R stars (Eriksson et al., Paper III), M stars (Plez et al., 
Paper IV), S stars (Plez et al., Paper V), C stars (J0rgensen et 
al., Paper VI) and models for very metal-poor stars (Paper VII). 
In these subsequent papers, the most important new opacity data 
are detailed, the model structures are described and analysed, 
and the model SED:s, fluxes and colours are compared to some 
observational data in order to explore and illustrate the applica- 
bility of the models. 

The total number of models is about 10 4 . The details of 
these models, including their calculated fluxes at different wave- 
lengths, are available via http://marcs.astw.uu.se 

2. The development of MARCS 

Since the early 1970ies, we have developed and used the 
code MARCS for constructing late-type model atmospheres. 
Spectral-line blanketing was first considered using Opacity 
Distribution Functions (Gustafsson et al. 1975) and this tech- 
nique and its underlying so-called ODF approximation was 
proven to be fully adequate for F, G and K stars. Thus, mod- 
els could be calculated with just a few hundred frequency points 
and extensive grids for G and K giants (Bell et al. 1976) as 
well as R stars (Olander 1981) were issued on the basis of the 
ODF approximation. The ODFs were constructed with a pro- 
gram called SSG, originally devised for synthetic colour calcu- 
lations by R.A. Bell (see, e.g. Bell 1971). The line list used by 
this program was based on laboratory wavelengths and oscillator 
strenghs, as well as "astrophysical oscillator strenghts" derived 
from the solar spectrum. The total list contained about 50,000 
lines. The SSG program was next used to calculate model spec- 
tra and colours, which were systematically compared with obser- 
vations (e.g. Gustafsson & Bell 1979; Bell & Gustafsson 1989). 
From these comparisons a generally good agreement was found, 
which made us confident in applying the models in a number 
of studies, e.g. of abundances in globular cluster stars. However, 
there was a mismatch in certain wavelength regions, in particular 
in the ultraviolet and violet, where the models were obviously 
too bright. We tentatively ascribed this discrepancy to the ab- 
sence in the line lists and the models of many weak lines, which 
together affect the stellar spectra. 

There were also more fundamental problems with the ODF 
method. One was that the ODF:s had to be recalculated when 
the chemical composition was changed, which made it inflexi- 
ble for stars with special and peculiar abundances, such as car- 
bon stars. To circumvent this, a special method to add ODF:s 
for individual molecules was invented (Saxner & Gustafsson, 
1984). However, for N-type carbon star models we also found 
that the effects of molecules in the upper atmospheric layers, 
contributing strong opacity at wavelengths different from the 
atomic absorption at greater dephts, made the ODF approxima- 
tion unsatisfactory (cf. Ekberg et al. 1986). To provide satisfac- 
tory analysis for a study of the chemical composition of car- 
bon stars (Lambert et al. 1986), we consequently changed the 
scheme to Opacity Sampling, in which the adequate absorption 
at each monochromatic wavelength point is treated in full detail 
(yet, assuming LTE in our case), following Peytremann (1974) 
and Sneden, Johnson & Krupp (1976). In order to obtain an ac- 
curate representation of the total radiation contribution to the 
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heat balance and the radiative pressure force this required an 
increase of the number of frequency points by one to two or- 
ders of magnitude. A number of model grids were calculated 
with this new program, for carbon stars (cf. Lambert et al. 1986; 
J0rgensen, Johnsson & Nordlund 1992), for M giants (Plez, Brett 
& Nordlund 1992) and M dwarfs (Brett & Plez, 1993) and ap- 
plied to studies of individual stars. We gradually extended the 
underlying line list by calculating and adopting more complete 
data for diatomic molecules, such as TiO (J0rgensen 1994a, 
Plez, Brett & Nordlund, 1992, Plez 1998), as well as polyatomic 
molecules like HCN (J0rgensen et al, 1985, J0rgensen 1990), 
C2H2 (J0rgensen thesis, cited in J0rgensen 1989) and H2O (Plez, 
Brett & Nordlund 1992, Alvarez and & Plez 1998, J0rgensen et 
al. 2001, Decinetal. 2000). 

In order to provide models for a major effort to study the 
build-up of chemical elements in the Galactic Disk using solar- 
type stars, we further extended the line list, by now adding the 
very great number of metal lines from Kurucz's calculations 
(Plez, Brett & Nordlund 1992; Edvardsson et al. 1993). This ver- 
sion of the code was also used in calculating a grid for hydrogen- 
poor carbon stars (R Cr B stars, Asplund et al. 1997), and ap- 
plied in an abundance study of such stars (Asplund et al. 2000b). 
Somewhat different versions of the code were used to compare 
models to ISO spectra and other IR spectra, partly for calibration 
purposes of the satellite instrumentation, of F, G and K stars (van 
der Bliek, Gustafsson & Eriksson, 1996, and subsequent papers, 
Decin et al. 2003 and references therein) as well as of M stars 
(Alvarez, Lanon, Plez et al. 2000. Fluks et al. 1994) and carbon 
stars (Loidl, Lancon & J0rgensen 2001). We also made com- 
parisons of model colours with observed ones (Bessell, Castelli 
& Plez, 1998), and applied colours and SEDs for establishing 
temperature scales (see, e.g., Massey et al. 2007 and references 
therein). 

The development of MARCS has always been driven by our 
own (and collaborators') needs for models for studies of par- 
ticular stars, mainly for abundance determination. It has been 
a flexible "laboratory set-up", rather than a "common-user in- 
strument" , and we have been reluctant to calculate models for 
stars of types that we cannot explore ourselves to understand 
the limitations of our models. Until now, less than 10 3 models 
have been published, which is a small number compared to the 
number actually calculated and used. Also, several different ver- 
sions of MARCS have appeared, which has led to some confu- 
sion because of inhomogeneities among the grids produced. It 
was therefore judged important to develop a common updated 
version, and to use it to construct and publish an extensive ho- 
mogeneous grid. The results of this effort is presented in this and 
subsequent papers. 



3. Physical assumptions and equations 

The basic assumptions have been listed above. We shall not 
discuss their adequacy further here (for a review, see e.g. 
Gustafsson & J0rgensen 1994), but comment on them one at a 
time, and in connection with that list some of the corresponding 
fundamental equations, for further reference below. 



3.1. Stratification and hydrostatic equilibrium 

Assuming spherical symmetry we may write the equation of hy- 
drostatic equilbrium: 



GM, 

VFtot = -p—r 



(i) 



Plot being the total pressure, p the matter density, G Newton's 
constant of gravity and M r the stellar mass inside radius r. We 
neglect the atmospheric mass in comparison with the total stellar 
mass M and thus assume M, = M. For VP to t we have 



V/\ot = VF„ + VAurb + VP ra 



(2) 



where P g and f tur b are the gas pressure and turbulent pressure, 
respectively, and the last term is the force excerted by the radia- 
tion, 



1 r 

c Jo 



(k a +cr A )F A dA, 



(3) 



where F A is the radiative energy flux per wavelength unit, k a and 
o- A are the monochromatic absorption and scattering coefficients, 
respectively, and c is the speed of light. The boundary condition 
for Eq. © is 



PJr = oo) = P tulb (r = 00) = 0.0 . 



(4) 



The formulation of the boundary condition involving the radia- 
tive force is, however, not quite trivial in practice; see Plez, Brett 
& Nordlund (1992). For the turbulent pressure f t uib one may 
write 



Pturb =jSpv t 



(5) 



where p is the gas density and v t a characteristic velocity. This 
pressure is measuring the force produced by the kinetic move- 
ments of the gas, whether due to convective or other turbulent 
gas motions. The parameter/? is — 1, with an exact value de- 
pending on whether the motions occur more or less isotropically. 
With our general lack of knowledge about v t it is reasonble to 
adopt an approximate recipe for P tur b- We begin with assum- 
ing a depth independent value of v t . The dominant depth vari- 
ation in P g = UpT ln mo \, % being the gas constant, occurs in p. 
Neglecting the variation in the temperature T and in the mean 
molecular weight yU mo i as well as in v t and in M, = M, one finds 



GM 



1 + «v? 



(6) 



where xf is the flux-weighted mean of the extinction coefficient 
per gram and L the luminosity of the star. Here, we have also 
assumed that the dominating fraction of the stellar flux is carried 
by radiation, i.e. the approximation is most accurate for the up- 
per radiative zones of the models. One may regard the right-hand 
side of Eq. (O as the local effective surface gravity g e ff times p. 
We thus find that we can mimic the effects of the radiative force 
and/or the turbulent pressure on the models by using models with 
those effects neglected with an adjusted gravity: 



eff - . 



where 



1 v? 



GM 



(1) 



(8) 



Defining the effective Eddington luminosity, L^ dd = 4ttGMc/xf 
we obtain 



( 1 - (L/Lf M ) ^ 



rfeff - . 



i+/?r(v t /c s ) 2 J 



(9) 
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where c s = ^yP g jp is the sound spead for an ideal gas and y 
is the adiabatic index. Basically, we have neglected the depth 
variation in y(v t /c s ) 2 in deriving this expression. 

Thus, a model with a turbulent velocity v t may be repre- 
sented by a model with a reduced gravity g e g and v t = 0, ac- 
cording to this recipe. Similarly, effects from the radiation force 
may be mimicked by changes in g or M. However, the Planck 
mean may vary strongly with depth and with stellar fundamen- 
tal parameters, and one should therefore be careful in the use 
of Eq. (0 for exploring the radiative effects on the atmospheric 
structure and, e.g., the possible effects on mass-loss rates. For 
a more detailed study for red supergiants and AGB stars, see 
Gustafsson & Plez (1992) and J0rgensen & Johnson (1992), re- 
spectively. We have tested the use of Eq. (|7]l to simulate the ef- 
fects of turbulent pressure for a number of models at various 
points in the grid and find that it leads to very small errors in the 
temperature structure (less than 5 K in the temperature through 
all the model for a depth independent v t in the interval to 10 
km/s. We therefore have chosen to set v t = for all grid mod- 
els, and advice those who would have liked a different choice to 
use models with a different mass or g, according to the recipe 
given in Eq. (0. It should be noted that the mixing-length treat- 
ment of convection adopted here (see Sec. 2.2 below) leads to a 
rapidly varying formal convective velocity, in particular close to 
the boundary where Schwarzschild stability sets in. If this varia- 
tion were included the term VP t urb = V(/3 p ) would get a ma- 
jor contribution from Vv t . However, more realistic simulations 
of convection show that v t varies much less with depth than p, in 
accordance with what was assumed in the derivation of Eq. (O 
above. 

The acceleration of gravity, g, for our spherical models is 
a depth-varying quantity, according to Eq. ([8j. Also the stellar 
energy flux F(r) and r e fj are varying with depth. Thus, we have 



where 



F(r) = F md (r) + F conv (r) = cr SB T eff (r) 



(10) 



where F rdd and F conv are the radiative and convective flux, re- 
spectively, and cr SB is Stefan-Boltzmann's constant. What re- 
mains constant (in stationary models) is the mass M (since only a 
tiny fraction of the total stellar mass resides in the atmosphere), 
and the luminosity, L = Anr 2 F(r). We label the models by the 
values of T e s and g at a radius r — R\ where tr oss , the optical- 
depth scale based on the Rosseland mean opacity, is equal to 1 .0, 
i.e.: 



T e ff = T e fi(R\) - 



(11) 



i + yipxRoss ef 

8 cr SB T 3 pXRoss C 



d\nT 

Vt = , 

dlnP 

and Hp is the local pressure scale height, 



H 



P 

gp' 



(15) 



(16) 



(17) 



V a d is the adiabatic temperature gradient, C p the specific heat at 
constant pressure, xros* the Rosseland mean opacity and I the 
mixing length. Eq. (fT~3T > is valid if (Vt - V a d) > 0. Also 



i GM 

-a/-^p2<Wv 



with 



* P \8T 



(18) 



(19) 



the derivative taken at constant thermodynamic pressure. In ad- 
dition to the mixing length parameter, a = £/H p , there are two 
more explicit parameters in this formulation: y which is related 
to the adopted temperature distribution within the convective el- 
ements and v which deals with the energy dissipation by the 
turbulent viscosity. As a standard, we have chosen the param- 
eters according to the suggestions by Henyey et al.: a - 1.5, 
y = 0.076 and v = 8. It should be noted that different choices 
than these are made, sometimes without being pointed out, in 
other current work. The effects of varying the convective param- 
eters were explored for models of red giants by Gustafsson et 
al. (1975). However, these variations do not at all map the real 
possible range of errors due to our (inadequate) treatment of the 
convective energy transport and even less can give relevant infor- 
mation on the effects of the thermal inhomogeneities generated 
by convection. 

The convective flux is added to the radiation energy flux, and 
we may then write the energy equilibrium equation 



+ F t 



ad 



4nr 2 ' 



(20) 



3.3. LTE and radiative transfer 



MG 



(12) 



3.2. Mixing length convection 

We have used the version of the Mixing-length "theory" as pre- 
sented by Henyey et al. (1965). The convective energy flux is 
given by 



All number densities of all atoms and molecules are assumed 
to follow from the corresponding laws for thermal equilibrium, 
the Saha equation and the corresponding equation of chemical 
equilibrium for the molecules (Gibson & Heitler 1928, see also 
Russell 1934). Similarly, all excitation equilibria and all parti- 
tion functions are calculated adopting the Boltzmann distribu- 
tion, with higher terms in the atomic partition functions cut ac- 
cording to the method of Irwin (1981). 

The radiation source function is assumed to be 



i e 

^CODV ~ TP T ^conv ~7T~ ■ 

I Mr, 



(13) S A = 



KA -B A (T) + 



KA + CTA 



-Ja, 



(21) 



Here, 



6 A : 



(l+D 



(V T - V ad ) , 



(14) 



where B^{T) is the Planck function. All line absorption is as- 
sumed to occur in true absorption, i.e. 



Ka 



cent + Hne 



(22) 
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where K° ont and a -1 ™ 6 are the sums of all continuous absorption 
contributions and line absorption contributions, respectively. 
The mean intensity Jx is calculated from 



J» 



Jo 



Qx)dfi. 



Here, following Feautrier (1964), we define jx(p) as 



(23) 



(24) 



where Ix(p) is the specific intensity in the direction specified by 
fj. = cosfl, 6 being the angle relative to a stellar radius. The equa- 
tion of radiative transfer for j,\ is 



d 2 h 



(25) 



where tx is measured along the ray. For a discussion of ade- 
quate boundary conditions for jx and their implementation, see 
Nordlund (1984). Knowing jxip), the wavelength-integrated flux 
in the radial direction can be calculated from a derivative of the 
second Eddington moment Kx, defined by 



Jo 

Thus, we have for the monochromatic flux 

\ OTx r k a + crx ) 
cf. Mihalas (1978, his Eq. 2.80) and then 

F m d = Jf a cIA, 

which is used in the energy conservation equation, Eq. (1201 1, 
Alternatively, the energy balance may be expressed as 



(26) 



(27) 



(28) 



—((F ad + F com )-r 2 ) = 0. 



(29) 



For F conv = which is often the case in the upper layers of the 
models, this is equivalent to 



J k a [Jx(tr oss ) - Bx(T(t Ro s S ))] dA = 

— ?rad — ^thermal — , 

with 



(30) 



'rad = J Kx Jx < 
nal = J Kx B A i 



'/rad = J Kx Jx(tr oss ) dA 

Q thermal = I KxBx(T(T Ross ))dA . 



(31) 



4. Physical data 

The volume of physical data needed in the calculation of model 
atmospheres of late-type stars is considerable. Data are needed 
for the calculation of the ionization equilibrium of atoms and dis- 
sociation equlibrium of molecules. These data include chemical 
composition data, ionization energies and dissociation energies, 
as well as partition functions. Moreover, continuous absorption 
and scattering coefficients are needed. The by far most extensive, 
and most demanding need, however, is the various data that are 



necessary for the proper calculation of the line absorption. Here, 
we shall briefly present the data used to calculate ionization- 
dissociation equlibria as well as continuous absorption and scat- 
tering. Some more details concerning the line-absorption data 
are given in Papers II - VII. 

The basic chemical composition adopted is that of the Sun, 
as listed by Grevesse, Asplund & Sauval (2007). There is still 
some dispute on the C, N and O abundances adopted there 
(C=8.39, N=7.78 and 0=8.66) and therefore we have alterna- 
tively taken the data of Grevesse & Sauval (1998) with CNO 
abundances higher by about 0.2 dex for the solar-metallicity 
sub-grid and explore the differences in Sec. 6.3. In varying the 
overall metallicity of the models ([Me/H]) we have as a stan- 
dard scaled the abundances of all elements heavier than He in 
unison. The abundances of so-called a elements O, Ne, Mg, Si, 
S, Ar, Ca, and Ti tend to vary somewhat differently for galac- 
tic stars (see, e.g. Edvardsson et al. 1993, Reddy et al. 2003, 
Ryde & Lambert 2004, Nissen et al. 2004, Cayrel et al. 2004); 
a variation which, however, is not necessarily the same for all 
stellar populations in the Galaxy (see, e.g. Fuhrmann 1998 and 
Bensby et al. 2005) or in other galaxies (see Venn et al. 2004). 
Therefore, we offer models with two different sets of abun- 
dances: one with a uniform scaling for each element with [Me/H] 
and one where the a elements have been scaled as follows: 
[or/Fe]=0.4 for -5.0 <[Me/H]< -1.0, [a/Fe]=-0.4x[Me/H] for 
-1.0 <[Me/H]< 0.0, and [a/Fe]=0.0 for [Me/H]> 0.0. For the 
giant stars, the results of the First Dredge Up of CNO processed 
material is known to lead to a reduced carbon abundance and a 
correspondingly increased N abundance (cf. e.g., Boothroyd & 
Sackmann 1999). We therefore also offer models with revised 
C and N such that C/N=1.5 and 0.5 by number, respectively, 
as alternatives to the solar value of C/N=4.07, though keeping 
C+N constant (cf Sec 6.3). For the corresponding giant stars, the 
changed C/N ratio is usually accompanied by a decrease in the 
isotopic ratio 12 C/ 13 C. For the "CN processed" models we have 
therefore changed this ratio from solar to 20 and 4, respectively. 
Finally, for the M stars and carbon stars we also present models 
with increased C abundances, such that C/O ranges from 0.54 to 
0.99 (M and S star models) to above 1.0 (models for C stars). 

The atomic partition functions and ionisation energies used 
are from Irwin (1981) with some modifications. The molecular 
equilibria were calculated using the partition functions and dis- 
sociation energies of Sauval (private communication), which is 
an update of Sauval & Tatum (1984) for diatomic molecules, 
and Irwin (1988, and private communication) for polyatomic 
molecules. Altogether 519 molecules were included in the equi- 
librium calculations. A comparison with equlibria calculated by 
Piskunov (private communication) shows a very good agree- 
ment. 

The continuous absorption from Hi, H , H 2 , H£, Hei, 
He",Ci, Cn, C", Ni, Nn, NT, Oi, On, O , Mgi, Mgn, All, 
Aln, Si i, Sin, Cai, Can, Fei, and Fen, as well as CH, OH, 
CO and H2O were calculated from sources according to 
Table 1 . Corrections were made to the data for C i, Mg i, Al i, 
Si i in order to account for the fact that TOP base photon 
cross sections are published with theoretical energy levels de- 
rived from the model atom, which depart significantly from 
the much better known data from laboratory work. For these 
four species corrections were considered necessary since they 
significantly affect the UV fluxes of the models. We identi- 
fied the individual TOP base energy levels for these species 
and shifted them to the empirical energies given by NIST 
( h ttp : //physics. n ist. go v/PhysRefDa ta/A SD/index. html). The TOP 
base gives the photon absorption cross section (in Mbarn) as 
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Table 1. Sources of data for continuous opacities, "b-f" and 
"f-f" denotes bound-free and free-free processes, respectively, 
and "CIA" is Collision Induced Absorption (or quasi-molecular 
absorption). "TOPbase" refers to the TOP data base of data 
from the Opacity Project (Seaton et al. 1994), available via 
http://vizier.u-strasbg.fr/topbase/home. html 



Ion and process 


Reference 


Hi b-f and f-f 


Karzas & Latter (1961)* 


H b-f 


Wishart (1979) 


H f-f 


Bell &Berrington (1987) 


He i b-f 


TOPbase 


He i f-f 


Peach (1970) 


He" f-f 


lohn (1994) 


Ci, Cnb-f 


TOPbase 1 


C i, C ii f-f 


Peach (1970) 


c-f-f 


Bell etal. (1988) 


Ni,Nnb-f 


TOPbase 


N-f-f 


Ramsbottom et al. (1992) 


Oi, On b-f 


TOPbase 


Of-f 


John (1975a, 1975b) 


Mg i, Mg ii b-f 


TOPbase f 


Mgif-f 


Peach (1970) 


Ah, Alnb-f 


TOPbase* 


Si i Si ii b-f 


TOPbase f 


Si i f-f 


Peach (1970) 


Cai, Can b-f 


TOPbase 


Fei b-f 


TOPbase, Bautista (1997) 


Fe ii b-f 


TOPbase 


other metals f-f 


Peach (1970), hydr. approx. 




Mihalas (1964) 


H- f-f 


John (1975a, b), John & Williams (1975) 


CH b-f 


Kurucz etal. (1987) 


OH b-f 


Kurucz etal. (1987) 


co- f-f 


John (1975a, 1975b) 


H z O- f-f 


John (1975a, 1975b) 


Hi+HiCIA 


Doyle (1968) 


Hi+HeiCIA 


Gustafsson & Frommhold (2001) 


H 2 +HiCIA 


Gustafsson & Frommhold (2003) 


H 2 +H 2 CIA 


Borysow et al. (2001) 


H 2 +HeiCIA 


J0rgensen et al. (2000) 


H 1 scattering 


Dalgarno, quoted by Gingerich (1964) 


H 2 scattering 


Dalgarno & Williams (1962) 


e~ scattering 


Mihalas (1978) 


He 1 scattering 


Bues&Wehrse(1976) 



' modified according to the occupation probability formalism, see text 

1 http://cdsweb. u-strasbg.fr/OP. htx 

The calculated wavelengths of UV and blue C 1, Mg 1, 

Al 1 and Si 1 absorption edges have been shifted to their 

laboratory values, see text 



a function of the energy of the ejected electron. This elec- 
tron energy was converted to the energy of the absorbed pho- 
ton by adding the ionization energy from the specific elec- 
tronic level. The cross-section data for all levels was summed 
up assuming LTE for a number of temperatures and plotted 
at the full original wavelength resolution at different temper- 
atures. Also collision-induced absorption of H1+H1, Hi+Hei, 
H2+H1, H2+H2 and H2+Hei was included. H2-He and H2- 
H2 CIA data from Borysow et al., referred to in Table 1, 
are available from http:www.stella.nbi.dk/pub/scan. More ex- 
tensive data sets for the CIA are publicly accessible from 
http: //www. astro, ku. dk/~aborysow . 

Continuous electron scattering as well as Rayleigh scatter- 
ing from H 1, H2 and He 1 was included and was assumed to be 
isotropic. 



Line absorption data for atoms and first ions were origi- 
nally adopted from VALD-1 (Piskunov et al. 1995) but recently 
data for the ten most important species were modified accord- 
ing to VALD-2 (Stempels, Piskunov & Barklem, 2001). The 
updated species are Sii, Cai, Can, Tin, Cri, Cm, Mnn, Fei, 
Fen, and Nin. The VALD-data adopted include the very nu- 
merous line data calculated by Kurucz & Bell (1995, see also 
http://kurucz.harvard.edu), including lines between predicted en- 
ergy levels that are not yet experimentally verified. The gf 
values were, however, modified for 948 lines from VALD-1 
and 1276 lines from VALD-2 in the wavelength range 3783 to 
8968 A, based on fits of synthetic solar intensity spectra. These 
model spectra were calculated with the Holweger & Miiller 
(1974) model and compared with the observed solar disk-center 
FTS spectrum of Brault & Neckel (1987), see Paper II. The 
model spectrum was required to match the observed equiva- 
lent widths to within 0. 1 dex in log gf .The resulting gf values 
are listed at http://m arcs.astro. uu. se\ Line absorption data for di- 
atomic and polyatomic molecules were considered for species 
according to Table 2. Most line lists of J0rgensen and collabo- 
rators are described by J0rgensen (1997) and may be obtained 
by anonymous ftp via www.stella.nbi.dk/pub/scan. Line lists for 
HCN, C2H2 and C3 will be added later to this data base in con- 
nection with the paper on carbon-enriched stars in the present 
series. The line lists of Plez referred to in Table 2 as "unpub- 
lished" are available via h ttp ://ma r cs.astro.uu.se\ Some of these 
lists were especially tailored for the present grid; in such cases a 
more detailed description of the line list for the particular species 
is to be found in one of the subsequent papers in the present se- 
ries. For the OH A-X system we did not use the data of Gillis 
et al. (2001) since the Kurucz (1995b) list includes higher vibra- 
tion and rotation states. For MgH we did not include the new 
data of Skory et al. (2003) for 24 MgH since we were anxious to 
be consistent with the MgH lines for other Mg isotopes. 

Hydrogen line and bound-free opacity, and their 
merging, were modelled using a code by Barklem 
(http://www.astro.uu.se/~barklem/hlinop.html), which is based 
on the occupation probability formalism of Dappen, Anderson 
& Mihalas (1987). The details of the description of this line 
opacity are given by Barklem & Piskunov (2003). For the 
atomic lines of metals, the damping wings were also calculated, 
using the best available data (Anstee & O'Mara, 1995, Barklem, 
Piskunov & O'Mara, 2000a and references therein, 2000b, 
Barklem & Aspelund-Johansson 2005, and Barklem, private 
communication). For weak lines where such data were missing 
we just adopted the Unsold recipe (cf. Unsold, 1955, his Eq. 
76.43), with a constant enhancement factor set to 2.0 for Nai, 
1 .3 for Si 1, 1.8 for Ca 1, 1 .4 for Fe 1 and 2.5 for all other species. 
(The damping constants given in the line list of Kurucz, which 
were obtained by perturbation theory, might have been preferred 
to the use of Unsold values.) For the molecules, damping was 
not taken into account - the lack of proper damping parameters 
makes any such attempt questionable. Also, the huge density 
of molecular lines often makes the damping wings of a line 
less important since the intensity of the Doppler cores of 
neighbouring lines dominate. 

In general, it should be noted that although the present 
atomic and molecular data are very extensive, and much more 
complete and accurate than, e.g., a decade ago, they are far 
from perfect. Considerable improvements are still needed, both 
by adding many more faint lines and determining line strengths 
more accurately. The risk that present models are still underblan- 
keted must be appreciated. 
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Table 2. Sources of data for molecular line opacities. 



Species 


Reference 


HCN vib-rot 


Harris et al. (2002), 




J0rgensen et al.(2001) 


H 2 vib-rot 


Barber et al. (2006) 


C2 Phillips, Swan, 




Ballik-Ramsay 


Querci et al. (1971, priv. comm.) 


c 3 


J0rgensen et al. (1989) 


C 2 H 2 


J0rgensen (1989) 


CH vib-rot 


J0rgensen et al. (1996) 


CH A-X, B-X, C-X 


Plez et al. (2008) 


CN A-X, B-X 


Plez (unpublished) 


CO vib-rot 


Goorvitch (1994) 


CO A-X 


Kurucz (1995) 


CaH A-X, B-X 


Plez (unpublished) 


FeH F 4 A-X 4 A 


Plez (unpublished) 


MgH A-X, B/-X 


Kurucz (1995a) 


NH A-X 


Kurucz (1995a) 


OH vib-rot 


Goldman et al. (1998) 


OH A-X 


Kurucz (1995a) 


SiH A-X 


Kurucz (1995a) 


SiO vib-rot 


Langhoff & Bauschlicher (1993) 


TiO a,/3, y, yi, 6, e, <p, 




E 3 n-B 3 n, a'A-f'A 


Plez (1998) 


VO A-X, B-X, C-X 


Plez (unpublished) 


ZrOB'n-A'A, B'n-X'I, 





C'l-X'EE'O-A'A, 
b 3 n-a 3 A, d 3 0>-a 3 A, 

e 3 n-a 3 A, f 3 A-a 3 A, Plez et al. (2003) 



5. Numerical methods 

5.1. The general method 

Equations ([T][2][3][l3l|20]or[30l|27]and|28l» with relevant bound- 
ary conditions form a closed system of equations for determin- 
ing the dependent variables F conv , T, P g , P m d, Ptmb, Pe, as well 
as quantities like v conv , F ra( j, or alternatively q m i. In addition to 
these quantities, we also need to calculate quantities characteris- 
tic (in LTE at least) of the local temperature and pressure: k(X), 
<t(A), H p , Q, V a( j and y. The system of equations is solved on 
a tross scale, with a-r oss calculated by integration over all ~ 10 s 
wavelength points. This choice of depth scale leads to temper- 
ature structures that are only moderately affected by changes in 
the fundamental stellar parameters. Basically, a standard multi- 
dimensional Newton-Raphson method is applied to solve the 
highly non-linear system. All equations are first discretized, in 
the variables tr oss and A. Next, all the resulting equations are 
linearized in the dependent variables listed above. The input 
variables at the calculation of thermodynamic quantities and ab- 
sorption coefficients are T and P e , which is a very appropriate 
choice as long as H is a dominating opacity source. (For the 
coolest models, e.g. Ph or P g could have been more advanta- 
geous.) As a result of the linearization, a set of linear equations 
in the variations of the dependent variables is obtained. The co- 
efficients of this system contain derivatives of a great number 
of quantities relative to the dependent variables. First, a starting 
model is assumed, in order to enable a first calculation of the 
coefficients. Next, the system of linear equations is solved nu- 
merically, the resulting variations are applied to the dependent 
variables, new coefficients are calculated and a new solution is 
obtained for the variations. The elimination scheme is basically 
that of Rybicki (1971), also applied by Gustafsson & Nissen 
(1972) and Gustafsson et al. (1975), with an elimination work 
that scales linearly with the number of wavelength points, n^. 




Fig. 1. Sketch of the geometric representation used for solving 
the spherically symmetric radiative transfer problem. The full 
radius mesh used for representing the variables as functions of 
radius (about 7 points per decade in r r ) is indicated along the 
radius vector by a set of concentric circles. The distance of rays 
from the centre of symmetry is measured by the impact param- 
eter p. About 15 rays with impact parameters < R c , the core 
radius, are used and are here represented by two dashed lines. 
The set of rays that do not hit the core, typically about 50 rays 
but dependent of wavelength, are here represented by five par- 
allel rays. The mesh of points used for solving the equation of 
radiative transfer are the crossing points between the circles and 
the horizontal lines in the figure. Distances along the rays are 
measured by z — pr where p -cos 9. 

5.2. Radiative transfer 

The total radiative flux, F ra< i, q rd & and V,f, a( j are non-local func- 
tional of the model structure and can be calculated from a so- 
lution of the transfer equation, Eq. (1251 1. For the model structure 
obtained after each iteration, we need an accurate solution of 
the equation with the source function S x given by Eq. (l2T1 > and 
Eq. d23l >. which then enables us to calculate F m d and q ra $ for 
that structure. These quantities are used to calculate the correc- 
tions (right hand sides) in the next Newton-Raphson iteration. 
We also need a linearization of F la( j and q m a in T and P e to 
calculate factors, symbolically written as e.g. 5F la( i(T,)/5r(T^), 
or 8q m i(ji) 1 6P e (jk), to be used in the coefficient matrix in the 
forth-coming iteration. Here, r, and are any two radial optical 
depths in the stellar atmosphere. I.e., we need to estimate how 
the radiation contribution to the energy balance in each point in 
the atmosphere is affected by variations in temperature or pres- 
sure anywhere else. Both these tasks are accomplished with the 
iterative method of Nordlund (1984) which was applied to model 
atmosphere calculations by Plez, Brett & Nordlund (1992) for M 
giants and by J0rgensen, Johnson & Nordlund (1992) for carbon 
stars. Here, we shall only give a brief summary for reference. 

The equation of radiative transfer is solved along a set of 
parallel rays, which cover a suitable range of impact parameters, 
chosen individually for each wavelength (cf. Fig.[T|i. In practice, 
about 15 rays equidistant in p (the cosine of the angle relative 
to the normal) are chosen for impact parameters < R c , where R c 
is the radius for the stellar non-transparent "core". This core is 

defined by a radial J > 200, in order to allow a great 

enough thermalization depth for scattered light. Rays that never 
hit this core but only pass through the transparent "atmospheric" 
region, are all chosen to go through pre-selected tr oss points in 
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the vertical scale along a line through the stellar centre and per- 
pendicular to the ray. These rays are distributed with about 5 rays 
per decade in Tj along this line. 

Scattering is handled by an iterative technique where the ra- 
diative transfer equation along a particular single ray is used to 
correct residual errors. This technique uses a "form factor" fx, 

Mr) = Ja/Ja- (32) 

Here, both J\ and jfo can be taken from the previous iteration, 
since their ratio is very little dependent on the detailed shape of 
the source function; in this respect the method is analogous to 
the variable Eddington technique of Auer & Mihalas (1970) in 
handling the angular dependence of the radiative field. After a 
few iterations a correct solution is obtained. When linearizing 
the energy equilibrium equation (Eq.[20lor[30l) we linearize the 
transfer equation along such a representative ray and apply the 
form factor, as calculated in the preceeding iteration. This relates 
the changes in the energy balance to changes of the specific in- 
tensity along the ray. We adopt this relation as typical for all rays. 
The final result is then a coupling of the energy equilibrium at 
each point in the atmosphere to temperatures and electron pres- 
sures all over the model. 

The representative ray is chosen in the following way: At the 
radius in the discretization of the Tr oss scale where ta measured 
radially is about ^{kj + cta)/ka the set of rays hits the sphere at 
a set of angles relative to the normal. We take the ray where the 
cosine for this angle is close to 1 / V3 - a choice which is natu- 
ral in view of a result of the Eddington approximation that the 
outgoing intensity from a gray atmosphere with a linear source 
function in this direction has a characteristic depth of formation 
which can be taken as an average depth of formation for all disk 
rays. 

In order to obtain a sufficient numerical accuracy we lin- 
earize the radiative flux directly (Eq.l20b only for optical depths 
T Ross > 0.01, while for smaller depths we linearize the diver- 
gence of the flux (i.e. Eq.[30b. This avoids the loss of accuracy as- 
sociated with the nearly constant monochromatic fluxes at most 
wavelengths in the optically thin layers. 

In the discretization in tr oss we use 56 points, distributed 
between tr oss = 10~ 5 and tr oss = 100. 

The integrals over A are calculated using the trapezian rule. 

In the A discretization we use 10,000 points during the first 
iterations. In the final iterations there are 108,000 points, set to 
produce a sampling density R = A/5A = 20, 000, with A m \ n = 
9 10 A for models warmer than 3500 K and 1300 A for cooler 
models, and A max = 20 pm. We have made a series of numerical 
experiments to investigate the errors resulting from the sampling 
of the spectrum with a resolution smaller than what corresponds 
to the characteristic line widths. This investigation has been 
made by making a number of sets of models with 30,000 wave- 
length points which are identical except for the set of wavelength 
points used. We then find the temperature errors for tr oss > 10 
to be less than 0.003xr e ff K with the maximum errors in the 
upper layers of the model (above tr oss = 0.01). The correspond- 
ing errors in the flux integrals and flux derivative integrals (i.e. 
the integrals in Eqs.d28l) and ( f30l > by the sampling procedure are 
less than 1.5%. For the grid models, with 108,000 wavelength 
points, the sampling errors are correspondingly smaller (by ap- 
proximately a factor of (108, 000/30, 000)" 1/2 ). They then lead 
to temperature errors of about 0.002xr e ff K or less. For the M 
and C star models, we have found the corresponding errors in 
fluxes and temperatures to be reduced by a factor of three or 
more. This is a result of the overlapping molecular-line absorp- 
tion which reduces the dominance of a few sampling-demanding 



strong metal lines in the flux derivative integrals in the upper 
model layers. 

It should be noted that the fluxes, even if they are monochro- 
matically precise within the model constraints, may be far from 
representative for the regions surrounding each wavelength point 
in the spectrum. Thus, model fluxes integrated to represent nar- 
row spectral regions, AA wide, will have relative errors of typ- 
ically 0.3 x [A/(RAA)] l/2 . E.g., for lOOAwide bands at 5000A, 
the sampling errors will correspond to standard deviations of 
typically 1.6%. Obviously, synthetic spectra more detailed than 
the model fluxes must be calculated for theoretical calibra- 
tions of e.g. narrow-band photometry. We have included the 
surface fluxes in all 108,000 wavelength points into the model 
files published on http://marcs.astro.uu.se - more detailed high- 
resolution synthetic spectra will also be gradually added. 

For the integration of the p moments of the radiative intensity 
(in practice ja) we use the p points defined by the angles in the 
crossing points between the different parallel rays with impact 
parameter /?,■ and the concentric spheres around the stellar centre 
with radii r k corresponding to the discretization of the vertical 
t Ross scale, T k cf. Fig 1. I.e., p ik = {r\ - pj) 1/2 /r k , directions for 
which jxip) is directly available. Considering the character of the 
radition field in a spherical atmosphere this is an adequate set of 
p points, in spite of the fact that they are different for different 
radii, and different wavelenghts. The typical number of p points 
is 6 in the plane-parallel case and typically 40-60 in the spherical 
models. The quadrature in p is performed using piece-wise cubic 
splines, as decribed by Nordlund (1984). 

5.3. Representation of line opacity 

The line opacities have been pretabulated in the following way: 
For each individual species (e.g. Hi, Ci, Cn, CH, C2 etc.) and 
the full set of wavelength points, arrays of line absorption were 
calculated for a number of T - P(, points (or T points for the 
molecules). Pg is given by 

P 6 = P(H I) + 0.42P(He I) + 0.85P(H 2 ) (33) 

where the static dipole polarizabilities and mean velocities of 
He 1 and H2 relative to that of H 1 have been used to approxi- 
mately scale the broadening for Hi to these species. Altogether, 
17 T values and 18 values were chosen. Such tables were 
contructed for 6 different values of the microturbulence param- 
eter, £ t = 0, 1, 2, 5, 10 and 20 km/s. Next, for a given chemi- 
cal composition and microturbulence parameter the tables were 
summed for all atoms and ions, over the different species, with 
due consideration to the ionization equlibria. This gave one ta- 
ble of the summed atomic line absorption at each lambda point 
for each T, Pf, and xi t chosen. For each molecular species a ta- 
ble was produced for a combination of 12 T and 6 £ t values. 
In the subsequent model -atmosphere calculations, the logarith- 
mic line absorption coefficients were found from these tables 
by spline interpolation to the appropriate temperatures and log- 
arithmic pressures. If models were needed at microturbulence 
parameters different from those of the tables, also interpolation 
can made. 

In the calculation of opacity data tables the relevant absorp- 
tion for each spectral line at each T and P(, value was added to 
the OS table for wavelengths points progressing to the blue and 
the red side of the central line wavelength as far as the line opac- 
ity exceeds a pre-defined cross section. Also for very weak lines 
each line was then included at at least two wavelength points. 
To limit the total computing time needed and still guarantee that 
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no significant opacity be lost, these limiting cross sections were 
empirically determined for groups of atoms and ions. 

5.4. Computing time, starting model and convergence 

The total computing time needed for the calculation of a model 
atmosphere (the pre-calculation of the absorption coefficient ta- 
bles excluded) scales approximately as follows: (1) For the ra- 
diative transfer part of the problem the time is proportional to 
rip x n T x nx- (2) For the calculation of ionization equilibria and 
molecular equilibria the time scales as n T and (3) for the calcula- 
tion of absorption coefficients as n T x n^. In addition to that, (4) 
some initiation time is needed. Setting N T = n T /28, Na = «,i/10 5 
and Nf, = n^/6 and noting that we usually set — 6 in the 
plane-parallel case and that w 2/3 x n T in the spherical case, 
we find empirically that the time f„ needed per iteration is, in 
seconds, 

t it ~ 10 + 1(W T + 5N A + 10N T N A + INrNxN^. (34) 

For a typical spherical model with 108,000 wavelength points 
about 80 seconds per iteration are needed on a Mac Pro quad 
Intel Xeon 2.66GHz computer for one processor. For a plane- 
parallel model this is reduced to about 60 seconds. The relative 
small reduction reflects the fact that interpolation in the line- 
absorption tables (the N t Na term) constitutes a major fraction 
of the computing time; only for a highly spherical case (when 
N M becomes great) the radiative transfer calculations dominate. 
No doubt, the calculations of absorption coefficients and their 
derivatives could be speeded up further, e.g. by pre-tabulation. 
This would primarily be of interest at the calculation of exten- 
sive model grids, or of models with more complex physics than 
decribed by the approximations in Sec. 3, e.g. as regards hydro- 
dynamics. 

In multi-dimensional Newton-Raphson schemes like the 
present one, convergence is rapid, provided that a starting solu- 
tion which is close enough to the final solution has been chosen. 

When model calculations are started from scratch with the 
MARCS program, usually a gray starting model is chosen for 
the radiative zone. When the calculated flux in the convec- 
tively unstable zone of the first starting model exceeds the to- 
tal flux, the temperature gradient is automatically and directly 
reduced, beginning at the onset of the convective instability, un- 
til the convective flux is smaller than the prescribed total flux. 
In the present grid, we usually start from a nearby model in the 
fundamental-parameter space. If the starting model has a differ- 
ent effective temperature, a simple scaling of the temperature 
structure r(TR OSS ) may be applied, but this is not necessary if 
steps of only a few hundred K are taken in T e ff. 

The route towards convergence is often rapid but not al- 
ways quadratical. The pronounced non-linearities, not the least 
in the temperature dependencies of the molecular equilibria with 
strong effects on the opacities may slow down the convergence 
if one is not very close to the final solution. Also, the chang- 
ing presence and depth of the convective zone with temperature 
and metallicity affect the convergence. Usually, convergence to 
temperature corrections althrough the models of less than a few 
K, and logarithmic pressure corrections smaller than 0.01 dex, is 
obtained after 4 to 10 iterations. For the models with r e jf > 4000 
K the convergence is close to quadratic with corrections reduced 
by almost one order of magnitude from one iteration to the next. 
For the cooler models the convergence is slower but still fast. 
For some parameter choices, however, converged models were 
not obtained. This is in particular the case for the models with 
a strong radiative pressure gradient (i.e. close to the effective 



Eddington limit, see Gustafsson & Plez 1992) which is close 
to upsetting the hydrostatic equilibrium condition, i.e. models 
with high temperatures and low gravities. Another difficulty oc- 
curs for models in a small band in the T e ff - logg diagram, ex- 
tending from about T eff = 8000/T, logg = 5, to T eS = 615QK, 
logg = 3, where convergence is not achieved in the deepest lay- 
ers of the models because the convectively unstable regions are 
swapping between thin convective zones and zones extending to 
depths below the bottom of the model. This difficulty may be cir- 
cumvented by extending the depth scale of the model to deeper 
layers. Some of the models, though not fully converged in the 
deepest layers, are nevertheless presented in the grid since the 
spectrum forming regions (above tr oss = 10) are not affected 
at all. Convergence problems occurred for some of the coolest 
models with H2 convection zones at the surface as well as for 
models where the radiative pressure totally dominates the gas 
pressure. 



6. General properties of the models 

All models of the grid are accessible via http://marcs.astro.uu.se/ 
The organization of the model files is described there. Details 
of structures, thermodynamic variables, molecular partial pres- 
sures and fluxes in 10 5 wavelength points are given in the files. 
In Fig.|2]we present some sample temperature structures of grid 
models with different effective temperatures and metallicities. It 
is seen that the variation with metallicity at the surface and in 
the deep layers increases when T e s decreases from 8000 K to 
5000 K, and that this variation with metallicity changes sign at 
4000 K in a non-trivial way. Subsequently, this and other prop- 
erties of the grid models will be commented on. Our discussion 
here is confined to general properties of the model structures, 
while more details on models of particular types of stars, as well 
as discussion of model fluxes are deferred to later papers in this 
series. 

6.1. General effects of blanketing 

The effects of spectral lines on stellar atmospheres ("blanket- 
ing") have been explored thoroughly since it became possible 
to calculate models with lines included in the 1960s and 1970s. 
Here, we shall make a short summary of these results as a back- 
ground for some more specific comments on the blanketing ef- 
fects in the new grid. 

In general, the blocking of radiation leads to heating of the 
deeper layers of the atmospheres ("back-warming"). Here and 
subsequently we measure the back-warming effect at a given 
continuum optical depth, in practice at A = 500 nm (i.e. not 
tr oss which for the cooler models is strongly affected by line ab- 
sorption). The back-warming is in fact easily estimated, directly 
from the definition of the effective temperature, by assuming that 
a model with effective temperature is heated such that the 
deeper layers correspond to an unblanketed model with an effec- 
tive temperature, 

T' eS = (l-X)-?.T eff , (35) 

where X is the fraction of the integrated continuous flux blocked 
out by spectral lines, 

poo 

X = ^_ . (36) 
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Fig. 2. The temperature structures for a set of model atmo- 
spheres with different T e g, logg = 3 and different metallicities. 



The validity of this relation has been tested for models with dif- 
ferent fundamental parameters by calculating X of these mod- 
els, and next unblanketed models with effective temperatures 
revised according to the relation were constructed. Next, these 
models were then compared with the corresponding blanketed 
ones. The test came out favourably - for models all over the 
parameter space of the grid, the r TRms structures in the interval 
- 1 < trqss < 0.5 were reproduced astonishingly well by the cor- 
responding unblanketed ones; the additional steepening of the 
temperature gradient accross the tr oss interval for the blanketed 
model with solar abundances may typically amount to 0.02 • Tem- 
per decade in tr oss . 

We have also calculated the total blocking fraction X for 
the grid models and find the remarkable behaviour shown in 
Fig. [3] Thus, X stays nearly constant for the temperature interval 
8000 K > r eff > 4000 K for a given metallicity, and is closely 
identical for dwarfs and giants, while below 4000 K it increases 
as temperature goes down, due to the heavy TiO and H2O block- 
ing. Here, the metallicity sensitivity of X decreseas (which re- 
flects the fact that the molecular lines fill in the continuum re- 
gions in the spectrum and define the tr oss scale; 7\tr oss ) being 
rather robust against metallicity changes) while the gravity sen- 
sitivity increases. 

For the upper layers of the atmospheres the situation is more 
complex, as was discussed e.g. by Gustafsson et al. (1975), 
Gustafsson & Olander (1979) and Carbon (1979) and as is il- 
lustrated in Fig. [4] 
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Fig. 3. The blocking fraction X in percent for models in the 
grid with two different metallicities. The dwarf models all have 
logg = 4.5 while the giant models have logg values increas- 
ing with temperature, from \ogg - 0.0 at T e g = 3000 K to 
logg = 3.0atr eff = 5000K. 



In short, for LTE the spectral-line absorption may cause a 
cooling or a heating, depending on whether the absorption is 
located on the long-wavelength side of the local Planck func- 
tion (cooling, Case 1) or on the short-wavelength side. In the lat- 
ter case, heating (2a) or cooling (2b) may result, depending on 
whether the absorption is concentrated to the upper layers (2a) 
or whether it extends to the deeper atmosphere (2b). Examples 
of Case 1 are the cooling by the IR VR lines of CO for models 
of late G stars (cf. Johnson 1973, Gustafsson et al. 1975) and 
by H2O in the IR for later M stars. A clear example of Case 
2a is the heating by the electronic bands of TiO in early M-star 
models (Krupp, Collins and Johnson 1978) while a Case 2b ex- 
ample is provided by the metal-line absorption in the UV and 
blue for F and G stars. We note in passing that these surface ef- 
fects are all dependent on the LTE assumption - in particular if 
the spectral lines are partially formed in scattering processes the 
effects may be significantly smaller due to the weaker coupling 
of the radiation to the temperature of the local gas. The back- 
warming effects, on the other hand, are not strongly dependent 
on the assumed line-formation mechanism as long as the amount 
of blocking is unchanged. Missing or erroneous line absorption 
data may, however, still be of significance as a source of system- 
atic errors in the back-warming effects. 

The general effects of blanketing are illustrated for the grid 
models by a test where all line absorption was swichted off. All 
models with T e g > 4000 K and solar metallicity show a similar 
reponse: the temperature is increased by about 10 % in the out- 
ermost layers if the line opacity is neglected, an increase which 
gradually diminishes to zero at tr oss = 0.01 and then changes to 
a lowering of the temperatures at depth since back-warming is 
inhibited, (cf. Fig.|4]i. For the coolest giants the effects become 
considerably greater. 

We shall now give some more specific comments on effects 
of blanketing in the new grid, with remarks on the effects of 
various opacity sources, and on the effects of changing abun- 
dances and microturbulence. Furthermore, we shall comment on 
the effects of sphericity, and explore some interesting coupling 
between these latter effects and blanketing. 
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Fig. 4. The effects of switching offline absorption on the temper- 
ature structure of a sequence of models with logg = 3.0 and solar 
metallicity. Note that AT = 7\nolines) - 7\lines). It is seen that 
the blanketing effects are fairly independent of effective temper- 
ature for models with r e ff > 4000. 

6.2. Effects of different opacity sources 

In Fig.[5]we present the effects of different opacity sources, mea- 
sured as integrated effects on the surface and at depth in a set of 
models for main-sequence stars. We have thus defined the quan- 
tities 



6 S = | Ar(T 5 )JlogT 5 / | <ilogT 5 

J-4 J-4 

6 b = J AT(T 5 )d\ogT 5 I J d\ogT 5 . 

AT is the quantity 

A7Xt 5 ) = [r(T 5 ) - T'(T 5 )VT eB 



(37) 



(38) 



where T(ts) is the temperature structure measured at the con- 
tinuum optical depth at 5000 A and T'(ts) the corresponding 
quantity for a modified model atmosphere with the same funda- 
mental parameters but for which line opacity sources have been 
neglected one at a time or all simulataneously (such that the cor- 
responding model is unblanketed). The reason why t$ and not 
tross were chosen here is that the Rosseland mean is directly 
affected by the spectral lines, in particular for the cooler stars, 
while we here wish to separate these effects on the temperature 
structure from those on the r scale. 

We see from Fig. [5] that for a sequence of main-sequence 
models with solar metallicity the total blanketing effect, both at 
the surface and in the deeper layers, is rather temperature inde- 
pendent. Thus, the mean surface cooling amounts to about 7% 
of T e ff and the mean heating of the layers below T5 = 10~ 2 is 
about 11%, over most of the effective-temperature interval. I.e., 
the somewhat decreasing importance of metal lines when pro- 
ceeding towards the hot end of our effective-temperature inter- 
val is rather well compensated for by the increasing hydrogen- 
line blanketing, and the diminishing effects of metals for the 
cooler models are compensated for by the molecular blanketing. 
These compensation effects explain why relatively old differen- 
tial abundance analyses, made by primitive model atmospheres 
such as scaled solar models, are often found to agree rather well 
with more recent results based on models with much more com- 
plete atomic and molecular data. In Fig.|6]we display the corre- 
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Fig. 5. Structural effects on the surface (black triangles) and 
depth (open circles) of model atmospheres by different groups of 
spectral lines. The solid lines show surface-heating and bottom- 
cooling as defined in Eq. (l37l i when all spectral lines are removed 
in the calculation of solar-metallicity dwarf models (logg = 
4.5). The dotted, dashed, and dash-dotted lines show, respec- 
tively, the effects of removing only the metal lines, the molecu- 
lar lines, or the hydrogen lines. The structural effects of spectral 
lines are surprisingly similar for all solar-type stars. The removal 
of all line opacity makes the models quite isothermal in the op- 
tically thin surface layers. 



sponding effects when absorption by different molecular species 
is left out. It is seen that the only significant molecular effects at 
the surface for these models are the surface cooling due to CO 
for moderately cool models and to H2O for the coolest ones and 
the surface heating due to TiO, while in the deeper layers the sig- 
nificant back-warming is due to TiO and H2O. Other molecules, 
like CN, play only a very minor role in this respect. 



6.3. Effects of abundance changes 

The variation of certain abundances in the MARCS grid lead to 
effects on the models. Such effects were explored for the CNO 
abundances of G and K giant models by Gustafsson et al. (1975), 
for M and S stars by Plez et al. (2003) and for carbon stars 
by Lambert et al. (1986) and J0rgensen, Johnson & Nordlund 
(1992). The results given here are illustrative and complemen- 
tary to those of previous studies - more detailed discussion of 
abundance effects will be presented in later papers in this series. 

The recent considerable downward corrections of solar CNO 
abundances by Asplund and collaborators (Grevesse, Asplund & 
Sauval 2007 and references therein) are found to lead to very 
minor effects for most regions in the parameter space spanned 
by the present grid. Thus, for the temperatures at given tr oss , the 
effects are less than 1 % x T e s for all depths and luminosities for 
giants at solar metallicity. The most pronounced effects in the 
surface layers appear in models with T e ff < 3500 K and models 
with r e fj ~ 5000 K where the reduced cooling due to H2O and 
CO, respectively, is visible. 
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Fig. 6. Similar to Fig|5] The structural effects of molecular lines 
separated for different species to show the results of removing 
the line opacity of TiO, water, and CO, respectively. 



The CNO abundances in stars are modified through the first 
dredge-up along the subgiant-giant branch, mostly converting C 
to N (see, e.g., Boothroyd & Sackmann 1999, Charbonnel 1994). 
The effects of this on the model atmospheres are, however, small 
or moderate, as illustrated in Fig. [7] The two most important ef- 
fects are again the enhanced H2O surface cooling for the coolest 
models (since more O is available if C and thus CO is depleted), 
and a diminished CO cooling for the somewhat hotter models. 
The changes due to effects by CN are marginal. The accompa- 
nying changes of 12 C/ 13 C ratios are obviously only of minor sig- 
nificance. 

The most severe atmospheric effects of CNO abundance 
changes occur, however, when the carbon abundance is raised 
by dredge-up on the AGB and, as the carbon abundance ap- 
proaches or exceeds that of oxygen, the star becomes an S or 
a C star, respectively, with drastic changes in radiation fields and 
thus temperature-pressure structure. This transition has been ex- 
plored using MARCS models (for references, see above), and 
will be further discussed in Paper V and Paper VI in the current 
series. As an illustration here, however, we present a sequence 
of models, with C/O ratios ranging from 0.5 to 2.4 in Fig. [8] 

The most metal-poor stars known until now show severe 
departures from "standard" CNO abundances: they seem to be 
comparatively very rich in the CNO elements (cf. Christlieb 
et al. 2004, Frebel et al. 2005, Norris et al. 2007 and refer- 
ences therein). In general, there is a group of low-metallicity 
stars with considerably enriched carbon abundances, so called 
Carbon-Enriched Metal-Poor (CEMP) stars, (see, e.g., Aoki et 
al. 2007), which may show carbon enrichments ranging in the in- 
terval between a factor of 10 and 1000 relative to a solar C/Fe, as 
well as similarly considerable N enrichments and possibly also 
O enrichments. The effects on the atmospheric structures from 
these enrichments are illustrated in Fig.[9]for one representative 
CEMP giant model. 

It is seen that the enhanced C and N abundances heat the 
layers -2 < logrs < 0.5 considerably for these model param- 
eters (a result which was also found by Hill et al. 2000). Also 
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Fig. 7. The effects of dredge-up of CN-cycled material to the 
surfaces of Pop I giant stars as illustrated by model atmospheres 
with different C and N abundances. In the top panel carbon has 
been converted to nitrogen such that the C/N abundance ratio has 
been lowered from the solar ratio of 4.07 to 1.5, typical for the 
first dredge-up in low-mass stars. In the lower panel the ratio has 
been further lowered to 0.5 and there also the effects of changing 
the 12 C/ 13 C ratio are displayed. The stronger back-warming of 
the 3000 K models is due to the increased importance of water 
as more oxygen gets available when the CO abundance is de- 
creased by the reduction of carbon while the heating effect at 
the surface merely reflects the inadequacy of T5 for this heavy 
line-blocked model - for a plot of the temperature vs. tr oss , an 
increased cooling at the surface appears. In the 5000 K model 
the decreased CO cooling leads to a hotter surface. 



the deep layers get warmer, which is due to the increased elec- 
tron contributions from the enhanced carbon. This increases the 
H opacity which lowers the density at a given optical depth and 
therefore reduces the convective energy flux which increases the 
inner temperature gradient. These effects change if also O is en- 
hanced correspondingly. Then the CO cooling takes over in the 
surface layers, and much of the C2 and CN heating absorption 
vanishes, since most of the carbon is bound in the CO molecules. 

The effects of the observed increase of [a/Fe] for galactic, 
though seemingly not all, metal-poor stars can also be explored 
using the grid models. Although these changes are easily seen 
in spectra of the models, e.g. in the strengths of the Mg and Ca 
lines, the effects on the temperature-pressure structures are only 
minor. As is seen in Fig. [10] they are, however, of some signifi- 
cance for models with < 4500 K where the higher Ti abun- 
dance increases the surface heating as well as the backwarming 
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Fig. 10. The effects of an increase in the a-element abundances by +0.4 dexon the temperature structure of models of dwarfs and gi- 
ants with [Fe/H]=-1.0. The comparison model in each case is a corresponding model with no ff-element enhancement ([a-/Fe]=0.0). 
The increased electron pressure lowers the density and convection efficiency at a certain Rosseland optical depth which increases 
the temperature gradient in the deeper parts of a model. 



of the TiO absorption. Simultaneously, the electron contributions 
from Mg and Ca increase and raise the H opacity, which partly 
reduces the effects of the Ti abundance change. As discussed 
above for the CEMP star models, the higher electron pressures 
also reduce the densities at a given optical depth, which dimin- 
ishes the convective flux and thus admits a stronger temperature 
gradient in the deep atmosphere. 

It should finally be noted that systematic errors in the mod- 
els, due to missing opacities, departures from LTE or convec- 
tively generated inhomogeneities and errors in calculated con- 
vective fluxes, may well vary systematically with chemical com- 
position. Thus, the trends discussed in the present section may 
have to be revised when more realistic models are available. 

6.4. Blanketing effects of turbulence 

In the present grid, there are models calculated with different 
microturbulence parameters % t , but with otherwise identical in- 
put data. This makes it possible to systematically explore the de- 
pendence of the blanketing effects on microturbulence. Before 
discussing the results some considerations of what one could ex- 
pect will be presented. 

An increased microturbulence enhances the blocking of the 
spectrum by isolated saturated spectral lines, in proportion to 

So, if all lines were saturated one might estimate from 
Eq. (l35l) that the backwarming AT/T would increase by about 
1 1 A ■ Atj,/{;, ■ X. However, this is an upper limit since a consider- 
able fraction of the spectral lines (though not necessarily those 
that contribute most significantly to the blocking) will not be sat- 
urated - their equivalent widths will be uneffected by a changing 
£ t and thus the total blocking not changed. Also, for more metal- 
rich stars the overlapping of strong lines reduces the differen- 
tial effects of a £ ( increase, as does the fact that many of these 
lines have strong damping wings. We thus expect the effects of 
microturbulence changes on the backwarming to be small for 
the most metal-poor stars where even the strongest lines are un- 
saturated, then increase but finally level off for the most metal- 
rich and coolest models, and never reach the estimated upper 
limit. These expectations are verified by the models, as is illus- 



trated in Fig.QT| The differential back-warming caused by an 
increased microturbulence parameter, e.g. by a factor of 2, is of 
about the same magnitude as the result of doubling the metal- 
licity. This suggests that the total blocking contribution of lines 
on the flat part of the curve-of-growth (mainly sensitive to 
is of similar significance as that of weak spectral lines (propor- 
tional in strength to the metallicity); this picture is complicated 
by both the effects of lines on the damping-part of the curve-of- 
growth (with only a square-root dependence on abundance) and 
the strong abundance sensitivity of some diatomic molecules. 
As regards the differential blanketing effects in the upper atmo- 
sphere when the microturbulence parameter is changed, the situ- 
ation is not quite obvious. The first, heating, part of the radiation 
equilibrium integral in Eq. ( f30b will not change in itself if only 
the line-absorption profiles become broader. If the effects on J,\ 
are taken into consideration, the decrease of the intensity over 
larger wavelength regions by the widening of the lines will, as 
long as they are weak, be compensated for by the increase of 
in the line center. However, for saturated lines will be reduced 
over broader wavelength intervals when is increased, and thus 
the heating term will become smaller. The second, cooling, term 
will not change from a broadening of the line profiles. Thus, 
a net cooling is expected to occur. We see from Fig. [TT] that in 
practice increasing from 1 to 2 km/s leads to some moderate 
cooling (on the order of 0.5 % in T/T e ff) of the upper layers for 
all effective temperatures while further increase up to 10 km/s 
leads to further cooling (by up to 2%, for the coolest models), 
however, this balance shifts to heating in the very surface layers. 

Changes in the macroturbulent parameter (classically repre- 
senting velocity shifts on geometrical scales corresponding to 
at least one optical depth unit in the continuum) are not able 
to change the model structures through radiative field changes; 
concerning effects of changes in the turbulent pressure which 
may be related to both macro- and microturbulent motions, see 
Sec 3, above. 

Real stellar atmospheres are not believed to be well repre- 
sented by the concepts of micro- and macroturbulence. Velocity 
gradients on intermediate scales are of vital importance, as has 
been shown by convective simulations (Nordlund & Dravins, 
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Fig. 8. A sequence of MARCS spherical model atmospheres 
with r eff = 3000 K, logg=0.0, M = 1 M Q . The models have dif- 
ferent abundance ratios, C/O, marked in the figure, but else solar 
metal abundances. Points on the temperature-pressure relations 
with logTRoss = -5, -4,-3,-2,-1,1 and 2 are marked with circles, 
which are filled for the carbon-rich models (C/0> 1.00). Points 
where logr Ross = are marked with x. The figure demonstrates 
the strongly increasing pressures at given temperature and opti- 
cal depth as the C/O ratio increases, until C/O reaches 1 .0 when 
the situation is drastically reversed as the opacities increase and 
the pressures decrease again. The models with C/O < 0.9 show 
M-type model spectra, those in the interval 0.9 < 0.98 S-type 
spectra, while those with greater C/O show N-star spectra. 



1990, Asplund et al. 2000a), and global gradients caused by 
pulsations as demonstrated in Lederer et al. (2006) and refer- 
ences therein. As regards the effects on the radiative field, these 
gradients allow for more heating in spectral lines since the gas 
high up in the atmosphere may absorb hot radiation from lower 
layers that was not absorbed below; similarly, the gas may cool 
more efficiently through spectral lines since the shifts make the 
gas more transparent in overlying layers. However, even more 
important are the hydrodynamical effects due to gas expansion 
and compression on the energy balance. Altogether, this makes 
the representation of the energy balance in real stars by classical 
models quite problematic; also, there is no hope to gain further 
insight into the physics of these systems by fine-tuning the tur- 
bulence parameters. 

6.5. Effects of sphericity 

The properties of spherically symmetric model atmospheres for 
late-type stars have been discussed or commented upon in nu- 
merous studies. The area was first pioneered by Schmid-Burgk 
& Scholz (1975) and Schmid-Burgk, Scholz & Wehrse (1981 
and papers cited therein), and later explored by Plez (1990), 
Plez, Brett & Nordlund (1992), J0rgensen, Johnson & Nordlund 
(1992), Hauschildt et al. (1999) and Heiter & Eriksson (2006). 
The papers by Plez and collaborators, by J0rgensen and collab- 
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Fig. 9. The effects of non-solar CNO abundance ratios on 
extreme Pop II star models with r e ff=4500K, logg = 1.0, 
[Me/H]=-3.0. In the top panel the carbon and nitrogen abun- 
dances have been increased by factors of 10, 100, and 300 above 
the solar ratio to iron, oxygen and all other elements. The heat- 
ing in the middle layers, primarily due to C2, is strong, and even 
more so at depth which is due to the increased electron pres- 
sures and continuous opacities. In the lower panel also oxygen is 
enhanced. The increased electron pressures in the more CN(O)- 
rich models lowers the density and convection efficiency at a 
certain optical depth which increases the temperature gradient. 
The major surface cooling agent in the outer surface layers is 
CO. 



orators and by Heiter & Eriksson, all build on earlier versions 
of the MARCS code. These papers demonstrate that the effects 
of sphericity are in general important for the temperature struc- 
tures, causing relative temperature effects on the order of 1% or 
more, for extensions of the atmosphere normalized on the stel- 
lar radius greater than about 2%. The relative extensions scale 
roughly as g~ 1 ^ 2 ; the sphericity effects are thus mainly important 
for low-gravity models. These papers do however not system- 
atically explore the effects on the models and fluxes of chang- 
ing from plane-parallel to spherical symmetric geometry. Here, 
some further comments will be made on this topic on the basis 
of the new model grid. 

The extension zq of the stellar atmosphere, measured as the 
geometrical thickness of the atmosphere (for instance defined 
as the layers with -5.0 < logT Ross < 0.0) is a key factor in 
the discussions. It is easy to estimate the extension from the 
hydrostatic-equilibrium equation, the definition of optical depth 
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Fig. 11. 7\tr oss ) for models with different microturbulence pa- 
rameters minus 7\tr oss ) for a model with £ t = 2 km/s, for so- 
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Thick lines denote models with £ t = 10 km/s, thin lines with 
£ t = 1 km/s. 



and the gas law for an ideal gas. The relative thickness of the 
atmosphere from optical depth r up towards the surface at to is 
thus: 



Z0 
R^ 



* <T> ln(-) = 4.10- 7 *?- r « . 

<Mmoi >gRi(n+ 1) r M g -(t7+1) 



(39) 



Here, we have assumed the variation of the pressure in the at- 
mosphere to be dominated by the density variation - the temper- 
ature and the mean molecular weight yU mo i in the gas law were 
thus approximated by a constant characteristic temperature and 
molecular weight, < T > and < /i mo i > (here set = 1.26), respec- 
tively. The mass absorption coefficient (opacity per gram) K g is 
assumed to vary (only) with the density p, K g = k-p' 1 . R® and M e 
are the model radius and mass, respectively, in solar units. 

According to our analytical estimate of the relative extension 
in Eq. ( 1391 ), the quantity 



C 



2.5- 10 6 [— ] 



Rf- 



(77+1) 



(40) 



' eff 



should be approximately constant and equal to 1.0. It is seen in 
Fig.[T2]that the relative extension is fairly well accounted for by 
the approximate expression, and that the choice of r\ — 1 is a 
reasonable fit for an extensive part of the HR diagram, which is 
to be expected since the dominating H absorption per gram is 
roughly proportional to the pressure. Note that the relative exten- 
sion of an atmosphere is relatively independent of the magnitude 
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Fig. 12. The value of C as calculated from Eq. ( |40l with model 
values of zo//?i for a set of models with r e jf ranging from 3000 K 
to 8000 K (with bigger symbols the higher the effective tem- 
perature) and metallicities [Me/H]=0.0, except for two models 
with [Me/H]= —1.0 which are attached with thin lines to corre- 
sponding [Me/H]=0.0 models, rj was chosen = 1. The value of C 
when based on the value of the extension zo/R\ estimated from 
Eq. ( 1391 is 1.0 and shown by the horizontal dotted line. 



of the opacity, which cancels in the zero-order approximation. A 
depth-variation of k is of some significance, but only marginally 
so. The most important sphericity effects related to the opacity 
are more indirect, through the effects that a changed temperature 
structure may lead to, in particular when molecular absorption is 
switched on, as will be demonstrated below. 

We now turn to the sphericity effects on the temperature 
structure of the models. The temperature in radiative equilibrium 
is set by the radiative field, following Eq. (f30b . The basic effects 
due to sphericity at a point high up in a stellar atmosphere come 
from the general dilution of the radiative field, partially due to 
the absence of incoming radiation from angles fi < fio where //o 
gradually increases the further from the optically thick layers the 
point is situated. Basically, this leads to less radiative heating of 
the upper layers of the atmosphere and thus cooler temperatures 
there. One may prove that this reduction of the mean intensity for 
small extensions will scale with (zo/R\) 1 ^ 2 . For larger distances 
from the star, a simple and approximative estimate for the tem- 
perature effect, was presented by Gustafsson et al. (1975) (see 
also Bohm-Vitense 1972) who assumed that the radiative flux 
should be roughly proportional to radius r~ 2 and to temperature 
T A as one may estimate from Stefan-Boltzmann's law. This leads 
to the estimate 



AT/T * -0.5 x zo/Ri 



(41) 



The effects of sphericity on the temperature structures of the 
models are shown in Fig.[T3]and compared with the approximate 
expression in Eq. (f4TT >. The agreement is satisfactory, though cer- 
tainly not perfect, which is partially explained by the approxi- 
mations made in deriving it. As important are the effects of the 
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Fig. 13. The difference in the surface temperature (at log tr oss = 
-5) between spherical and plane-parallel models with T e ff rang- 
ing from 3000 K to 7000 K (with bigger symbols the higher 
log g) and metallicities [Me/H] =0.0, except for two models with 
[Me/H]= -1.0 which are attached with thin lines to correspond- 
ing [Me/H] =0.0 models (the lower symbols in the pairs repre- 
senting the more metal-poor models). Three models at T e ff = 
4000 K have also been calculated with TiO absorption neglected 
and are denoted by crosses - it is seen that the positions of these 
models depart very significantly from the corresponding stan- 
dard models. The relation according to Eq. (|4TT > is shown as a 
straight line. 
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Fig. 14. The difference between temperature structures of spher- 
ically symmetric and corresponding plane-parallel models with 
logg = 1.0, M = 1.0 Mq, [Me/H] =0.0, and r eff ranging from 
3000 K to 6000 K. The complex variation of the sphericity ef- 
fects in the upper layers when TiO sets in as a heating agent at 
T ef{ * 4000 K is illustrated. 



wavelength dependence of k^, which weights the mean intensity 
and the Planck function B \ in the radiation equlibrium inte- 
grals of Eq. ( f30b differently for different models, depending in 
particular on the molecular absorption, and resulting in the large 
scatter in Fig.[T3l 

Now turning to the detailed model results, the sphericity ef- 
fects on AT/T are shown as a function of Tr oss for a sequence 
of supergiant models with different effective temperatures in 
Fig. [14] Obviously, the effects shift with T e ff, First being rela- 
tively small for the hottest models, they considerably increase 
for T e ff ~ 5000 K, then diminish again to reach a minimum at 
about 3800 K and finally again increase strongly as T e ff goes 
towards 3000 K. This behaviour may be interpreted as a cou- 
pling between the sphericity effects and the effects of molecu- 
lar blanketing. The coupling may show up as a positive or neg- 
ative feed-back. Sphericity cools the upper layers, relative to 
a plane-parallel model, and the molecular formation increases. 
Around 5000 K this leads to a strong surface cooling by CO. 
For cooler models TiO forms, which conversely heats the up- 
per atmosphere. Finally, H2O takes over as a dominating opacity 
source, and that again cools the model surface. That this explains 
the behaviour can be tested by blocking out different molecular 
opacities from the model calculations. In Fig.[l5]we thus show 
a sequence of models at r e tj = 4000 K with different log g and 



with and without TiO. It is clear that the differential sphericity 
effects above Tr oss are changed significantly by the TiO absorp- 
tion. From this figure we also see that the sphericity effects in the 
temperature structure get significantly greater than 1 % for mod- 
els with logg < 1 .0. However, if turbulent pressures are allowed 
for this value may be higher, according to Eq.©. 

7. Comparison with other model grids 

We shall here make some comparison of MARCS grid models 
with the models of other existing grids, but confine this compar- 
ison to the temperature structures - some additional comments 
on calculated fluxes will be be made in subsequent papers of this 
series. 

The most extensive and widely used grid of mod- 
els for most of our parameter space are the plane- 
parallel models of R. Kurucz and collaborators (available on 
http://kurucz.harvard.edu). Among the different sets of models 
published there, it is particulary relevant to compare with those 
that are calculated by Castelli & Kurucz 2003) with a standard 
mixing-length theory (without "convective overshoot"), which 
also seem to produce more consistent model parameters when 
applied to real data (Castelli, Gratton & Kurucz, 1997). Some 
sample comparisons of model structures are shown in Figs. [16] 
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Fig. 15. The difference between temperature structures of spher- 
ically symmetric and corresponding plane-parallel models with 
r eff = 4000 K and different log g. M - 1.0 M Q and [Me/H]=0.0, 
and with the TiO absorption alternatively switched on and off 
(dashed lines, the gravities of these can be traced via the mod- 
els to which they adhere at depth). The complex variation of the 
sphericity effects in the upper layers due to TiO is illustrated. 

[T7]and[T8j where we have chosen Castelli & Kurucz ODFNEW 
models to compare with. It is seen that the agreement in the tem- 
perature structures is almost perfect for the models of solar-type 
dwarf stars of different metallicities. This is true also for models 
of early M dwarfs. For the models of giants and supergiants, the 
agreement is again good, although the MARCS models tend to 
be some 10-80 K cooler in the surface layers from tr oss < 10~ 2 . 
The ageement is also very satisfactory if pressures or densities 
are intercompared. In view of the fact that these two grids of 
models are made with two totally independent numerical meth- 
ods and computer codes, with independent choices of basic data 
(although Kurucz's extensive lists of atomic line transitions are 
key data underlying both grids), this overall agreement is both 
satisfactory and gratifying. 

The situation was somewhat less satisfactory for models 
from the NextGen grid of spherically symmetric models, calcu- 
lated with the PHOENIX code (Hauschildt et al. 1999, models 
obtained by private communication in 1996). In the inner photo- 
spheres the agreement is good but again a difference systemat- 
ically appears at the surface, starting around tr oss w 10~ 2 , with 
the MARCS models beeing cooler, but now by as much as typ- 
ically 250 K. (As is seen in Fig. [7] the agreement between the 
Castelli & Kurucz and the MARCS models is generally much 
better than with the NextGen models.) However, more recent 
LTE models (PHOENIX- ACES -2008) obtained by Hauschildt 



Fig. 16. MARCS model atmospheres (solid) for dwarf stars with 
logg = 4.5, [Me/H]=0.0 and corresponding models from the 
Castelli & Kurucz ODFNEW grid (dashed). The curves are la- 
beled with relevant values of T e g. 



with an updated version of the PHOENIX code (Hauschildt, pri- 
vate communication) and which are also based on Kurucz line 
data, agree very well with the MARCS models. The basic rea- 
son for this improvement in the fit seems to be due to changes in 
molecular opacities. Again, the present overall fit is very sat- 
isfactory, in view of the independence of the methods of the 
MARCS and PHOENIX teams. 

A number of comparisons with models made by previous 
versions of MARCS have also been performed. In Fig. [20] we 
thus compare with models of G and K giants (Bell et al. 1976) 
calculated with the first version of MARCS. It is seen that the 
present temperature structures agree quite well with the older 
ones in spite of the much more incomplete line absorption data 
of those models. (An exception from this is seen in the outer 
layers of the 4000 K model, where the TiO absorption, missing 
in the BGEN models, contributes a surface heating of the more 
recent model). The pressures in the present models are, how- 
ever, significantly lower which reflects their more heavy opac- 
ities. The good agreement in the temperature structures must 
thus at least partially be fortuitous, e.g. reflecting that the line 
absorption missing in the 1976 models is not heavily biassed to- 
wards the ultraviolet or the red spectral regions relative to the 
absorption that is in. In Fig.|20]we also show one example where 
a spherically symmetric model for a late-type M supergiant is 
compared with a previous spherically symmetric MARCS model 
of Plez et al. (1992). The considerable differences in the upper 
photosphere are due to the new and improved molecular data 
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Fig. 17. MARCS model atmospheres for giants and supergiants 
with [Me/H]= 0.0 (solid) and corresponding models from the 
recent Castelli & Kurucz grid (dashed). The curves are labeled 
with relevant values of T e ff and log g. 

used now and its more detailed representation, in particular for 
the H2O absorption. 

8. Interpolation between the models 

It is useful to be able to interpolate the model structures to a 
set of parameters T s s, log g, and [Me/H] different from those 
tabulated in the grid. Such an interpolation is not completely 
straightforward, as strong non-linearities appear in the behaviour 
of some of the thermodynamic variables when stellar param- 
eters are varied. One program was developed by T. Masseron 
(2006) and extensively tested by him on a previous grid of 
MARCS models encompassing the following range of param- 
eters: 3800 K< r eff <7000 K, 0.0 < log g < 5.0, -4.0 < 
[Me/H] < 0.0. It interpolates T, P g , P e , kr oss , and the geomet- 
rical depth as a function of tr oss . The interpolated model must 
lie inside a complete cube of existing models in the parameter 
space (r e jf, log g, [Me/H]), and the interpolation is optimized 
to account for non-linearities in the grid. With the actual grid 
parameter steps, maximum errors in the interpolated quantities 
should remain below 0.25% for T and a few % for P g and P e , 
whereas a-r oss and the column density (rhox) are more difficult to 
interpolate. The program, with a detailed user manual, is avail- 
able on the MARCS web site: http://marcs.astro.uu.se 

An interpolation routine, applied to structures, as well as 
fluxes and colours of models by Kurucz and collaborators, has 
also been described by Nendwich et al. (2004), see also Valenti 
& Fischer (2005). 
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Fig. 18. MARCS model atmospheres with [Me/H] =-2.0 (solid) 
and corresponding models from the Castelli & Kurucz 
ODFNEW grid (dashed). The curves are labeled with relevant 
values of T e g and log g. 



9. Conclusions 

It will certainly take additional time before model atmospheres 
for late-type stars are regularly constructed with physically more 
adequate assumptions than the standard ones - of spherical sym- 
metry, mixing-length convection and LTE - adopted here. Until 
then, standard models will be the common choice in e.g. abun- 
dance analyses. Here we have presented such a grid of standard 
models to the state of the art. It includes model atmospheres for 
stars of spectral types from F to M, as well as carbon stars and 
Pop II stars of different metallicities. The assumptions, physi- 
cal data used and numerical methods have been described, and 
a number of general properties of the models have been illus- 
trated. In particular, we have discussed the effects of blanket- 
ing and sphericity in some detail and tested several of the exist- 
ing approximate analytical estimates of these effects and found 
them to work quite well. Also, already existing grids, such as 
the grid of plane-parallel models by Castelli & Kurucz and the 
new PHOENIX grid, have been compared with and the agree- 
ment has been found to be excellent in the region of overlap- 
ping parameters. This is particularly gratifying in view of the 
fact that the different approaches, as regards numerical methods 
and codes and to a considerable degree also the selection of data, 
are independent. We may then basically regard the calculation of 
standard ID LTE model atmospheres as becoming routine, but 
important improvements in line absorption data are still needed, 
and certainly not easily achieved. This is also true for data for 
dust formation and dust opacities. 
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Fig. 19. Model atmospheres with solar abundances, from the 
MARCS, NextGen and PHOENIX- ACES -2008 grids and (for 
r eff = 5000 K) the ODFNEW Castelli & Kurucz grid, logg is 
3.0 (5000 K) and 0.0 (3000 K), respectively. The MARCS and 
NextGen models are spherically symmetric with M = 5M Q 
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curve every decade in tr oss is marked by a o sign, and tr oss = 0.0 
by a x. The general impression of a better agreement between 
MARCS and Castelli & Kurucz models than between MARCS 
and NextGen models is also shown if the temperature is plotted 
relative to optical depth. 
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Further presentations and discussions of particular proper- 
ties of these models will appear in subsequent papers in this se- 
ries. Hopefully, these will serve also as starting points for further 
studies of the adequacy and shortcomings of the standard mod- 
els, and as a stimulus for endeavours to set higher standards in 
the art of modelling stellar atmospheres in the future. 
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